;********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/heat_stress.ncl"
;********************************************
 begin

;set baseline
 ;-----process one of the following 12 models each time in turn----
 ;mo  = (/"ACCESS-CM2","AWI-CM-1-1-MR","BCC-CSM2-MR","CMCC-CM2-SR5","EC-Earth3","GFDL-ESM4", \
 ;              "IITM-ESM","KIOST-ESM","MIROC6","MIROC-ES2L","MPI-ESM1-2-HR","MRI-ESM2-0"/); 12 models

 mo = (/"ACCESS-CM2"/); take this model for example

 syr = 1990
 fyr = 2099
 csyr = sprinti("%0.4i",syr);
 cfyr = sprinti("%0.4i",fyr);
 nyr = fyr - syr + 1;

;--data paths (input bias-corrected CMIP6 data; output 3-hourly G data; output daytime mean G for plots)
 dirPath = "~/CMIP6/models-vcbc/"+mo+"-vcbc/"     ; input bias-corrected CMIP6 data (*not in the repository, due to large data size ~15 TB)
 dirPath2 = "~/CMIP6/data-vcbc/3hourly_G/"   ; output 3-hourly G data (*not in the repository, due to large data size ~5 TB)
 dirPath3 = "~/source_data/outdoor/"  ; ouput daytime mean G (*provied in the repository)

;--read dimension for 1 deg simulations
 dimPath = "~/source_data/domain/"
 landfrac = "sftlf_fx_360x180_1deg_gn.nc"
 area = "areacella_fx_360x180_1deg_gn.nc"

 fdim  = addfile(dimPath+landfrac, "r")
 fdim2 = addfile(dimPath+area, "r")
 lat := fdim->lat
 lon := fdim->lon
 nlat = dimsizes(lat);
 nlon = dimsizes(lon);
 garea = fdim2->areacella  ;m^2
 lfrac = fdim->sftlf/100; %
 wgt := garea*1e-6*lfrac ; sqrkm
 wgt@_FillValue = lfrac@_FillValue; make sure FillValue is set correctly for wgt
 copy_VarCoords(garea, wgt)
 printVarSummary(wgt)
 print("global total land area: "+sum(wgt)+" sq km")



;====Human body Energy Balance Model (HEMB)====
;==constant parameters==
 alpha=0.3  ;reflectance of human body depending on clothes
 rho=1.2 ;density of air (kg/m3)
 cp=1010 ;specific heat capacity of air (Joules/kg/°C)
 ;=baseline skin temperature
 Ts=36 ;skin temperature in hot conditions, beyond which core temperature will rise (see Supplementary Table 1)
 esat_skin = satvpr_water_bolton(Ts, (/0,1/))  ;saturation vapor pressure at skin temperature 
 ;=radiation terms
 e_skin = 0.97 ;skin emissivity (Brewster, M. Quinn (1992))
 fs = 0.8  ;the fraction of total skin area effectively involved in radiant heat exchange (Steadman 1979, Fiala et al. 1999)

 ;=metabolic heat
 ;mean human skin surface area 1.7 (1.6-1.8 m^2; Parsons 1999)
 ;total metabolic heat generation 300 W for moderate work (Kjellstrom 2013)
 MH = 300/1.7  ;outdoor moderate work (W m^-2)
 

 DATA = True; process the model
 if (DATA) then

  nstep = 8; steps per day; 3-hourly data
  nhr = 24/nstep

  Gday_mx_sun = new((/nyr,nlat,nlon/),"float");
  Gday_mx_diff = new((/nyr,nlat,nlon/),"float");
  Gday_mx_diff_L = new((/nyr,nlat,nlon/),"float");
  Gday_mx_diff_U = new((/nyr,nlat,nlon/),"float");
  Gday_mx_zero = new((/nyr,nlat,nlon/),"float");
  Gnight_mx = new((/nyr,nlat,nlon/),"float");
  TWday_mx = new((/nyr,nlat,nlon/),"float");

  Gday0d_sun = new((/nyr,nlat,nlon/),"float"); annumal number of days with daytime average G>0
  Gday0d_diff = new((/nyr,nlat,nlon/),"float");
  Gday0d_zero = new((/nyr,nlat,nlon/),"float");
  Gnight0d = new((/nyr,nlat,nlon/),"float"); annumal number of days with nighttime average G>0
  TWday35d = new((/nyr,nlat,nlon/),"float");
   
  ;array of concurrent climate variables for Gday_mx (daytime average)
  Ta_Gday_sun = new((/nyr,nstep,nlat,nlon/),"float")
  RH_Gday_sun = new((/nyr,nstep,nlat,nlon/),"float")
  Rs_Gday_sun = new((/nyr,nstep,nlat,nlon/),"float")
  Rd_Gday_sun = new((/nyr,nstep,nlat,nlon/),"float")
  Rg_Gday_sun = new((/nyr,nstep,nlat,nlon/),"float")
  Ld_Gday_sun = new((/nyr,nstep,nlat,nlon/),"float")
  Lg_Gday_sun = new((/nyr,nstep,nlat,nlon/),"float")
  U_Gday_sun = new((/nyr,nstep,nlat,nlon/),"float")
  P_Gday_sun = new((/nyr,nstep,nlat,nlon/),"float")
  Rsa_Gday_sun = new((/nyr,nstep,nlat,nlon/),"float"); absorbed Rs
  Z_Gday_sun = new((/nyr,nstep,nlat,nlon/),"float"); zenith angle
  G_Gday_sun = new((/nyr,nstep,nlat,nlon/),"float")
  Ta_Gday_diff = new((/nyr,nstep,nlat,nlon/),"float")
  RH_Gday_diff = new((/nyr,nstep,nlat,nlon/),"float")
  Rs_Gday_diff = new((/nyr,nstep,nlat,nlon/),"float")
  Rd_Gday_diff = new((/nyr,nstep,nlat,nlon/),"float")
  Rg_Gday_diff = new((/nyr,nstep,nlat,nlon/),"float")
  Ld_Gday_diff = new((/nyr,nstep,nlat,nlon/),"float")
  Lg_Gday_diff = new((/nyr,nstep,nlat,nlon/),"float")
  U_Gday_diff = new((/nyr,nstep,nlat,nlon/),"float")
  P_Gday_diff = new((/nyr,nstep,nlat,nlon/),"float")
  Z_Gday_diff = new((/nyr,nstep,nlat,nlon/),"float")
  G_Gday_diff = new((/nyr,nstep,nlat,nlon/),"float")
  Ta_Gday_zero = new((/nyr,nstep,nlat,nlon/),"float")
  RH_Gday_zero = new((/nyr,nstep,nlat,nlon/),"float")
  Rs_Gday_zero = new((/nyr,nstep,nlat,nlon/),"float")
  Rd_Gday_zero = new((/nyr,nstep,nlat,nlon/),"float")
  Rg_Gday_zero = new((/nyr,nstep,nlat,nlon/),"float")
  Ld_Gday_zero = new((/nyr,nstep,nlat,nlon/),"float")
  Lg_Gday_zero = new((/nyr,nstep,nlat,nlon/),"float")
  U_Gday_zero = new((/nyr,nstep,nlat,nlon/),"float")
  P_Gday_zero = new((/nyr,nstep,nlat,nlon/),"float")
  Z_Gday_zero = new((/nyr,nstep,nlat,nlon/),"float")
  G_Gday_zero = new((/nyr,nstep,nlat,nlon/),"float")

  ;array of simultanous climate data for TWday_mx
  Ta_TWday = new((/nyr,nstep,nlat,nlon/),"float")
  RH_TWday = new((/nyr,nstep,nlat,nlon/),"float")
  Rs_TWday = new((/nyr,nstep,nlat,nlon/),"float")
  Rd_TWday = new((/nyr,nstep,nlat,nlon/),"float")
  Rg_TWday = new((/nyr,nstep,nlat,nlon/),"float")
  Ld_TWday = new((/nyr,nstep,nlat,nlon/),"float")
  Lg_TWday = new((/nyr,nstep,nlat,nlon/),"float")
  U_TWday = new((/nyr,nstep,nlat,nlon/),"float")
  P_TWday = new((/nyr,nstep,nlat,nlon/),"float");
  Z_TWday = new((/nyr,nstep,nlat,nlon/),"float");
  TW_TWday = new((/nyr,nstep,nlat,nlon/),"float");

  do ii=0,nyr-1;
    myr = syr+ii; model year

    print("calculate model: "+mo+" for year "+myr)
    f1C = systemfunc("cd "+dirPath+"; ls tas_3hr_"+mo+"_*_r*_g*_"+myr+"0101-*1231-vcbc.nc")
    f1 = addfile(dirPath+f1C, "r") ;distinguish from "addfile"
    f2C = systemfunc("cd "+dirPath+"; ls hurs_3hr_"+mo+"_*_r*_g*_"+myr+"0101-*1231-vcbc.nc")
    f2 = addfile(dirPath+f2C, "r")
    f3C = systemfunc("cd "+dirPath+"; ls ps_3hr_"+mo+"_*_r*_g*_"+myr+"0101-*1231-vcbc.nc")
    f3 = addfile(dirPath+f3C, "r")
    f4C = systemfunc("cd "+dirPath+"; ls rsds_3hr_"+mo+"_*_r*_g*_"+myr+"0101-*1231-vcbc.nc")
    f4 = addfile(dirPath+f4C, "r")
    f5C = systemfunc("cd "+dirPath+"; ls rlds_3hr_"+mo+"_*_r*_g*_"+myr+"0101-*1231-vcbc.nc")
    f5 = addfile(dirPath+f5C, "r")
    f6C = systemfunc("cd "+dirPath+"; ls rsdsdiff_3hr_"+mo+"_*_r*_g*_"+myr+"0101-*1231-vcbc.nc")
    f6 = addfile(dirPath+f6C, "r")
    f7C = systemfunc("cd "+dirPath+"; ls rlus_3hr_"+mo+"_*_r*_g*_"+myr+"0101-*1231-vcbc.nc")
    f7 = addfile(dirPath+f7C, "r")
    f8C = systemfunc("cd "+dirPath+"; ls rsus_3hr_"+mo+"_*_r*_g*_"+myr+"0101-*1231-vcbc.nc")
    f8 = addfile(dirPath+f8C, "r")
    f9C = systemfunc("cd "+dirPath+"; ls sfcWind_*3hr_"+mo+"_*_r*_g*_"+myr+"0101-*1231-vcbc.nc")
    f9 = addfile(dirPath+f9C, "r")

    dimtime := f1->time
    ;atts = getvaratts(dimtime)
    calendar = dimtime@$"calendar"$ ;get attribute value
    print("calendar type: "+calendar)

    if ( calendar.eq."360_day") then
      print("process 360_day data, year "+myr)
      nday = 360
    else if (calendar.eq."noleap" .or. calendar.eq."365_day") then
      print("process non-leap data, year "+myr)
      nday = 365
    else if (isleapyear(myr)) then
      print("process leap year "+myr)
      nday = 366
    else
      print("process noleap year "+myr)
      nday = 365
    end if
    end if
    end if
    print(" ")


    ;read annual data sequence for each model year
    ta := f1->tas
    rh := f2->hurs
    pa := f3->ps
    Rs := f4->rsds; total solar radiation at surface
    Ld := f5->rlds; atmosphere downward longwave at surface
    Lg := f7->rlus; upward longwave from ground
    Rd := f6->rsdsdiff; Solar downward near-infrared + visible diffuse at surface
    Rg := f8->rsus; surface_upwelling_shortwave_flux_in_air (reflected)
    ws10 := f9->sfcWind; reference height=10m
    ;printVarSummary(ta)

    taC := ta-273.15 ; K to degC
    copy_VarCoords(ta, taC)
    rh := where(rh.gt.100, 100, rh); reset supersaturation to 100% (mainly in high latitudes)
    rh := where(rh.lt.0, 0, rh); ;reset negative values (very few, not affect G in hot conditions)
    ws10 := where(ws10.lt.0, 0, ws10)
    Rs := where(Rs.lt.0, 0, Rs)
    Rd := where(Rd.lt.0, 0, Rd)
    Rg := where(Rg.lt.0, 0, Rg)   
    copy_VarCoords(ta, Rs)
    copy_VarCoords(ta, Rd)
    copy_VarCoords(ta, Rg)

    ;=method 1: only consider the top of a 2m object
    wind2m := 0.75*ws10 ;10-m wind to 2m wind FAO method
    delete([/ws10/])
    wind2m = where(wind2m .gt. 0.1, wind2m, 0.1); set minimal wind to 0.1 for natural convection (<0.1m/s is stale) 

    ;=first calculate solar constant and sun zenith angle at each time step
    ; Method used in CAM5.0 for Diurnal Cycle and Earth Orbit (tech note page 169-171)
    ; original citation Berger, A. L., Long-term variations of daily insolation and quaternary climatic changes,
    ; J. Atmos. Sci., 35, 2362–2367, 1978

   fSI = dirPath+"/SI_"+mo+"_"+myr+"0101-"+myr+"1231-vcbc.nc"
   if (.not. fileexists(fSI)) then

    ;--input constants to calculate sun zenith angle
    pi = get_pi(1)
    eps=23.4441*pi/180 ;epsilon: obliquity (deg to radians)
    ec=0.016715 ;eccentricity factor
    omega=(102.7 + 180)*pi/180 ;longitude of the perihelion + 180°(deg to radians)
    lam0 = 2*((ec/2+ec^3/8)*(1+sqrt(1-ec^2))*sin(omega) - \
               ec^2/4*(1/2+sqrt(1-ec^2))*sin(2*omega) + \
               ec^3/8*(1/3+sqrt(1-ec^2))*sin(3*omega) )

    ;--cd calendar day for each 3-hourly time step (0~364 with decimals) in UTC/GMT 0
    ;--use the timestamp of rsds (middle of 3-hour period) to estimate sun angle
    cdtime := f4->time   ; time:units = "days since 1850-01-01 00:00:00"
    cd := cdtime - floor(cdtime(0))    ; cd restarts at the beginning of each year
    print("calendar days: "+cd(0:10)); make sure it starts from 0 or 0.125 each year

    lam := lam0 + 2*pi*(cd-80.5)/365
    la := lam + (2*ec - ec^3/4)*sin(lam - omega) + \
         5*ec^2/4*sin(2*(lam - omega)) + \
         13*ec^3/12*sin(3*(lam - omega))
    delta := asin(sin(eps) * sin(la))
    t = 4
    print("solar declination on calendar day "+cd(t)+" = "+delta(t)*180/pi)

    ;--solar constant corrected by distance factor
    r0 := (1-ec^2)/(1+ec*cos(la - omega))
    SI := tofloat(1367*r0^-2) ;
    print("mean solar constant corrected by distance factor: "+avg(SI))
    SI := conform_dims(dimsizes(ta), SI, 0) ;convert to 3d [time, lat, lon]

    copy_VarCoords(ta, SI)
    fn0 = addfile(fSI,"c")
    fn0->SI = SI
   else
    fn0 = addfile(fSI,"r")
    SI := fn0->SI
   end if

   fzenith = dirPath+"/zenith_"+mo+"_"+myr+"0101-"+myr+"1231-vcbc.nc"
   if (.not. fileexists(fzenith)) then

    cd3d  := conform_dims(dimsizes(ta), cd, 0) ;convert to 3d [time, lat, lon]
    delta3d := conform_dims(dimsizes(ta), delta, 0) ;convert to 3d [time, lat, lon]
    lat3d := conform_dims(dimsizes(ta), lat, 1) ;convert to 3d [time, lat, lon]
    lon3d := conform_dims(dimsizes(ta), lon, 2) ;convert to 3d [time, lat, lon]
    Ha := 2*pi*(cd3d+lon3d/360) ;hour angle of sun during the day
    ;printVarSummary(Ha)
    coszen := sin(lat3d*pi/180)*sin(delta3d) - cos(lat3d*pi/180)*cos(delta3d)*cos(Ha)
    coszen := tofloat(coszen)
    zenith := acos(coszen)*180/pi ;solar zenith angle
    altitude:=90-zenith

    print("Range of coszen: "+min(coszen)+"~"+max(coszen))
    print("Range of zenith: "+min(zenith)+"~"+max(zenith))
    print("Range of altitude: "+min(altitude)+"~"+max(altitude))
    delete([/cd3d,delta3d,lat3d,lon3d,Ha/])

    copy_VarCoords(ta, zenith)
    fn1 = addfile(fzenith,"c")
    fn1->zenith = zenith
   else
    fn1 = addfile(fzenith,"r")
    zenith := fn1->zenith
    altitude:=90-zenith
    pi = get_pi(1)
    coszen := cos(zenith*pi/180)
   end if

  ;==calculate 3-hourly G for sun/shade/dark scenarios
  fGsun = dirPath2+"/G_sun_3hr_"+mo+"_"+myr+"0101-"+myr+"1231-vcbc-outdoor.nc"
  fGdiff = dirPath2+"/G_diff_3hr_"+mo+"_"+myr+"0101-"+myr+"1231-vcbc-outdoor.nc"
  fGzero = dirPath2+"/G_zero_3hr_"+mo+"_"+myr+"0101-"+myr+"1231-vcbc-outdoor.nc"
  if (.not. fileexists(fGsun)) then

    ;=use Fiala convective heat transfer function 
    hc := 14.1*wind2m^0.5 ;forced convection function from Fig.3 in Fiala et al. (1999)
    hc := where(wind2m .le. 0.1, 3.3, hc); set natural convection hc=3.3 W m-2 K-1 for U<=0.1 m/s

    ;=saturation vapour pressure of the air
    esat_air  := satvpr_water_bolton(taC, (/0,1/)) ;degC to Pa
    e_a := (rh/100)*esat_air

    ;=longwave radiation input (from air and ground) and output (emitted by skin)
    RLin := e_skin*0.5*(Ld + Lg)
    RLout = e_skin*5.67*10^-8*(Ts+273.15)^4 ;

    ;=latent heat of vaporization as a function of temperature==
    ;should use skin and sweat temperature, not air temperature
    lambda = latent_heat_water(Ts, (/0, 0/), 1, False)
    gamma0 := (pa*cp)/(0.622*lambda)

    ;==scenario 1: exposure to full solar radiation (outdoors)

    Rd:= where(coszen .gt. 0.1, Rd, 0);  Rd@_FillValue);
    copy_VarCoords(ta, Rd)
    print("mean incoming diffuse radiation: "+avg(Rd))
    print("mean incoming shortwave radiation: "+avg(Rs))

    ;--sun angle correction for direct beam according to Steadman (1979)
    ;1) convert direct beam radiation incident on a horizontal surface to that received on a surface perpendicular to the beam.
    Rbp := tofloat((Rs-Rd)/coszen) ;
    Rbp := where(coszen .gt. 0.1, Rbp, 0) ; only consider the day time
    Rbp := where(Rbp .lt. 0.9*SI, Rbp, 0.9*SI); constrain with solar constant, 0.9 is average fraction of direct beam in clear sky
    print("Range of Rb received normally perpendicular to beam")
    printMinMax(Rbp, 0)
    ;2) human projected area factor "phi" (Steadman 1979b)
    phi := 0.386-0.0032*altitude
    phi := where(altitude .gt. 85, 0.114, phi)
    phi := where(altitude .lt. 10, 0.354, phi)
    ;3)total absorbed shortwave with corrrection for direct beam-Rb, diffuse radiation-Rd, and reflected shortwave-Rg
    Rba := phi*(1-alpha)*Rbp 
    Rsa := phi*(1-alpha)*Rbp + 0.5*(1-alpha)*(Rd+Rg)
    Rsa := where(coszen .gt. 0.1, Rsa, 0) ; only consider the day time

    print("Range of incoming Rs")
    printMinMax(Rs, 0)
    print("Range of absorbed Rs")
    printMinMax(Rsa, 0)

    ;=full sun angle correction
    Rn := Rsa + RLin - RLout
    ;=add an effective radiation area factor (Steadman 1979a)
    Rn = fs*Rn ; apply the effective radiant area factor fs=0.8 (Steadman 1979a) 

    ;=scenario 2: exposure only to unavoidable diffuse radiation (under shade)
    Rn2 := RLin - RLout + 0.5*(1-alpha)*(Rd+Rg)
    Rn2 = fs*Rn2 ; apply the effective radiant area factor fs=0.8 (Steadman 1979a)

    ;print("mean incoming diffuse radiation: "+avg(Rd))
    ;print("mean reflected shortwave radiation: "+avg(Rg))
    ;print("mean net radiation at skin surface-diffuse: "+avg(Rn2))


    ;=scenario 3: exposure only to longwave radiation (completely dark)
    Rn3 := RLin - RLout 
    Rn3 = fs*Rn3 ;  apply the effective radiant area factor fs=0.8 (Steadman 1979a)


    ;==calculate energy balance G under different radiation scenarios 
    ;---sensible and latent heat fluxes
    H  := hc*(Ts - taC)
    LE := hc/gamma0*(esat_skin - e_a)

    ;1) for acclimatized persons 
    LEmax = 500 ; [W m-2] or 1.25 litre h-1 according to ISO 7933:2023
    ;2) for unacclimatized persons 
    ;LEmax = 400 ; [W m-2] or 1 litre h-1 according to ISO 7933:2023
    LE := where(LE .gt. LEmax, LEmax, LE);


    ;=scenario 1: full solar radiation
    Gs := Rn - H - LE + MH  
    print("min & max of Gs") 
    printMinMax(Gs, 0) ; 
    copy_VarCoords(ta, Gs)

    ;=scenario 2: diffuse radiation only
    Gdiff := Rn2 - H - LE + MH 
    copy_VarCoords(ta, Gdiff)
    print("min & max of Gdiff")
    printMinMax(Gdiff, 0) ;

    ;=scenario 3: zero solar radiation (longwave only)
    Gzero := Rn3 - H - LE + MH
    copy_VarCoords(ta, Gzero)
    print("min & max of Gzero")
    printMinMax(Gzero, 0) ;

    ;=save 3 hourly output
    system("bash -c 'if [ -e "+fGsun+" ]; then rm -f "+fGsun+"; fi'")
    fn1 = addfile(fGsun,"c")
    fn1->G_sun = Gs

    system("bash -c 'if [ -e "+fGdiff+" ]; then rm -f "+fGdiff+"; fi'")
    fn2 = addfile(fGdiff,"c")
    fn2->G_diff = Gdiff

    system("bash -c 'if [ -e "+fGzero+" ]; then rm -f "+fGzero+"; fi'")
    fn3 = addfile(fGzero,"c")
    fn3->G_zero = Gzero
  else
    fn1 = addfile(fGsun,"r")
    Gs := fn1->G_sun
    fn2 = addfile(fGdiff,"r")
    Gdiff := fn2->G_diff
    fn3 = addfile(fGzero,"r")
    Gzero := fn3->G_zero

  end if

  ;--use +-50% of Rd+Rg as upper/lower range of solar radiation under shade conditions
  ;-50% corresponding to mean transmittance of 16% under tree canopies
    Rsa1 := 0.5*(1-alpha)*(Rd+Rg)/2 ;
    print("mean lower limit of Rsa: 50% of Rd+Rg: "+avg(Rsa1))
  ;+50% corresponding to mean transmittance of 50% under shading nets
    Rsa2 := 1.5*(1-alpha)*(Rd+Rg)/2 ;
    print("mean upper limit of Rsa: 150% of Rd+Rg: "+avg(Rsa2))

    Rn_diff1 := RLin - RLout + Rsa1; lower
    Rn_diff2 := RLin - RLout + Rsa2; upper
    Rn_diff1 := fs*Rn_diff1
    Rn_diff2 := fs*Rn_diff2

    ;--lower/upper range of Gday_diff under the shade scenario
    G_diff1 := Rn_diff1 - H - LE + MH
    G_diff2 := Rn_diff2 - H - LE + MH
    copy_VarCoords(ta, G_diff1)
    copy_VarCoords(ta, G_diff2)

    opt = True
    opt@nval_crit = 1 ; require >=1 non-NA values per day to calculate day/night avg (default 1)
    Gdiff1_day := where(coszen .gt. 0, G_diff1, ta@_FillValue)
    copy_VarCoords(ta, Gdiff1_day)
    Gday_diff_L := calculate_daily_values(Gdiff1_day, "avg", 0, opt) ; daytime mean
    Gdiff2_day := where(coszen .gt. 0, G_diff2, ta@_FillValue)
    copy_VarCoords(ta, Gdiff2_day)
    Gday_diff_U := calculate_daily_values(Gdiff2_day, "avg", 0, opt) ; daytime mean


    ;==daytime/nighttime mean G under each scenario (daytime determined by solar zenith angle)
    opt = True
    opt@nval_crit = 1 ; require >=1 non-NA values per day to calculate day/night avg (default 1)
    ;----sun scenario---
    fGday_sun = dirPath2+"/Gday_sun_"+mo+"_"+myr+"0101-"+myr+"1231-vcbc-outdoor.nc"
    if (.not. fileexists(fGday_sun)) then
      Gs_day := where(coszen .gt. 0, Gs, Gs@_FillValue); set nighttime values to NA
      copy_VarCoords(Gs, Gs_day)
      Gday_sun := calculate_daily_values(Gs_day, "avg", 0, opt) ; daytime mean G (sun)
      fn1 = addfile(fGday_sun,"c")
      fn1->Gday_sun = Gday_sun
    else
      fn1 = addfile(fGday_sun,"r")
      Gday_sun := fn1->Gday_sun
    end if

    ;---shade scenario---
    fGday_diff = dirPath2+"/Gday_diff_"+mo+"_"+myr+"0101-"+myr+"1231-vcbc-outdoor.nc"
    if (.not. fileexists(fGday_diff)) then
      Gdiff_day := where(coszen .gt. 0, Gdiff, Gs@_FillValue); daytime only
      copy_VarCoords(Gs, Gdiff_day)
      Gday_diff := calculate_daily_values(Gdiff_day, "avg", 0, opt) ; daytime mean G (shade)
      fn1 = addfile(fGday_diff,"c")
      fn1->Gday_diff = Gday_diff
    else
      fn1 = addfile(fGday_diff,"r")
      Gday_diff := fn1->Gday_diff
    end if

    ;----dark scenario----
    fGday_zero = dirPath2+"/Gday_zero_"+mo+"_"+myr+"0101-"+myr+"1231-vcbc-outdoor.nc"
    if (.not. fileexists(fGday_zero)) then
      Gzero_day := where(coszen .gt. 0, Gzero, Gzero@_FillValue); daytime only
      copy_VarCoords(Gs, Gzero_day)
      Gday_zero := calculate_daily_values(Gzero_day, "avg", 0, opt) ; daytime mean G (dark)
      fn1 = addfile(fGday_zero,"c")
      fn1->Gday_zero = Gday_zero
    else
      fn1 = addfile(fGday_zero,"r")
      Gday_zero := fn1->Gday_zero
    end if

    ;---night time---
    fGnight = dirPath2+"/Gnight_"+mo+"_"+myr+"0101-"+myr+"1231-vcbc-outdoor.nc"
    if (.not. fileexists(fGnight)) then
      Gzero_night := where(coszen .le. 0, Gzero, Gzero@_FillValue); night
      copy_VarCoords(Gs, Gzero_night)
      Gnight := calculate_daily_values(Gzero_night, "avg", 0, opt) ; night mean
      fn1 = addfile(fGnight,"c")
      fn1->Gnight = Gnight
    else
      fn1 = addfile(fGnight,"r")
      Gnight := fn1->Gnight
    end if

    ;---wetbulb temperature---
    fTW = dirPath2+"/TW_3hr_"+mo+"_"+myr+"0101-"+myr+"1231-vcbc.nc"
    if (.not. fileexists(fTW)) then
      tdC := dewtemp_trh(ta, rh) - 273.15  ; C
      TW3h  := wetbulb(pa/100, taC, tdC)
      copy_VarCoords(ta, TW3h)

      system("bash -c 'if [ -e "+fTW+" ]; then rm -f "+fTW+"; fi'")
      fn1 = addfile(fTW,"c")
      fn1->TW = TW3h
    else
      fn1 = addfile(fTW,"r")
      TW3h := fn1->TW
    end if

    fTWday = dirPath2+"/TWday_"+mo+"_"+myr+"0101-"+myr+"1231-vcbc.nc"
    if (.not. fileexists(fTWday)) then
      TW_day := where(coszen .gt. 0, TW3h, TW3h@_FillValue)
      copy_VarCoords(TW3h, TW_day)
      TWday := calculate_daily_values(TW_day, "avg", 0, opt) ; 
      fn1 = addfile(fTWday,"c")
      fn1->TWday = TWday
    else
      fn1 = addfile(fTWday,"r")
      TWday := fn1->TWday
    end if


    ;========find annual maximum of daily max/avg/daytime/nightime means
    ;important to set negative missing values to sort in decreasing order
    ;only Gday_sun/Gday_diff/Gnight have missing values due to daylight variation
    TWday@_FillValue = -1.e+20
    Gday_sun@_FillValue = -1.e+20
    Gday_diff@_FillValue = -1.e+20
    Gday_diff_L@_FillValue = -1.e+20
    Gday_diff_U@_FillValue = -1.e+20
    Gday_zero@_FillValue = -1.e+20
    Gnight@_FillValue = -1.e+20

    print(nday) ;some data has leap years
    Gsort_sun := dim_pqsort_n(Gday_sun, -2, 0) ; Gsort_* only keeps the sorting sequence of original data
    Gsort_diff := dim_pqsort_n(Gday_diff, -2, 0) ;
    Gsort_diff_L := dim_pqsort_n(Gday_diff_L, -2, 0) ;
    Gsort_diff_U := dim_pqsort_n(Gday_diff_U, -2, 0) ;
    Gsort_zero := dim_pqsort_n(Gday_zero, -2, 0) ;
    Gnight_index := dim_pqsort_n(Gnight, -2, 0) ;
    TWsort_index := dim_pqsort_n(TWday, -2, 0) ; sorting x in decreasing order and return the sequence index

    ii_00 = 0 ;annual maximum of daytime mean
    Gday_mx_s = Gday_sun(ii_00,:,:)
    Gday_mx_d = Gday_diff(ii_00,:,:)
    Gday_mx_z = Gday_zero(ii_00,:,:)

    ;=save the hottest day ranked by each metric
    Gday_mx_sun(ii,:,:) = Gday_mx_s
    Gday_mx_diff(ii,:,:) = Gday_mx_d
    Gday_mx_diff_L(ii,:,:) = Gday_diff_L(ii_00,:,:)
    Gday_mx_diff_U(ii,:,:) = Gday_diff_U(ii_00,:,:)
    Gday_mx_zero(ii,:,:) = Gday_mx_z
    Gnight_mx(ii,:,:) = Gnight(ii_00,:,:)
    TWday_mx(ii,:,:) = TWday(ii_00,:,:)


   CALC = True;  False
   if (CALC) then

    ;===construct new day indices for each diel time steps
    ;day index ddd starting from 0 for 3-hourly G, to match with day index of Gday
    ddd := ndtooned(conform_dims((/nday,nstep/),ispan(0,nday-1,1), 0))

    ;==find concurrent 3-hourly values of climate variables when annual maximum daytime G happens==
    do i=0,nlon-1
      do j=0,nlat-1

        Gday1d = Gday_mx_s(j,i) ;
        if(ismissing(Gday1d)) then; skip ocean grids
           continue
        end if

        ip_00 = Gsort_sun(ii_00,:,:); index of original day for the Gday [nlat, nlon]
        ip = ip_00(j,i) ;day index of Gday varies across grids
        iday= ind(ddd.eq.ip); indices of 24hr/8 time steps for Gday
        if(all(ismissing(iday))) then; skip some wrong indices
           continue
        end if
        Ta_Gday_sun(ii,:,j,i) = taC(iday,j,i)
        RH_Gday_sun(ii,:,j,i) = rh(iday,j,i)
        Rs_Gday_sun(ii,:,j,i) = Rs(iday,j,i)
        Rd_Gday_sun(ii,:,j,i) = Rd(iday,j,i)
        Rg_Gday_sun(ii,:,j,i) = Rg(iday,j,i)
        Ld_Gday_sun(ii,:,j,i) = Ld(iday,j,i)
        Lg_Gday_sun(ii,:,j,i) = Lg(iday,j,i)
        U_Gday_sun(ii,:,j,i)  = wind2m(iday,j,i)
        P_Gday_sun(ii,:,j,i)  = pa(iday,j,i)
        Z_Gday_sun(ii,:,j,i)  = zenith(iday,j,i)
        G_Gday_sun(ii,:,j,i) = Gs(iday,j,i)

        ;=diffuse scenario
        Gday1d = Gday_mx_d(j,i) ;
        if(ismissing(Gday1d)) then; skip ocean grids
           continue
        end if
        ;= use permutation index of Gday_d
        ip_00 = Gsort_diff(ii_00,:,:); index of original day for the 99th percentile [nlat, nlon]
        ip = ip_00(j,i) ;day index 
        iday= ind(ddd.eq.ip); indices of 24 time steps for Gday
        if(all(ismissing(iday))) then; skip some wrong indices
           continue
        end if
        Ta_Gday_diff(ii,:,j,i) = taC(iday,j,i)
        RH_Gday_diff(ii,:,j,i) = rh(iday,j,i)
        Rs_Gday_diff(ii,:,j,i) = Rs(iday,j,i)
        Rd_Gday_diff(ii,:,j,i) = Rd(iday,j,i)
        Rg_Gday_diff(ii,:,j,i) = Rg(iday,j,i)
        Ld_Gday_diff(ii,:,j,i) = Ld(iday,j,i)
        Lg_Gday_diff(ii,:,j,i) = Lg(iday,j,i)
        U_Gday_diff(ii,:,j,i)  = wind2m(iday,j,i)
        P_Gday_diff(ii,:,j,i)  = pa(iday,j,i)
        Z_Gday_diff(ii,:,j,i)  = zenith(iday,j,i)
        G_Gday_diff(ii,:,j,i) = Gdiff(iday,j,i)


        ;=zero solar scenario
        Gday1d = Gday_mx_z(j,i) ;
        if(ismissing(Gday1d)) then; skip ocean grids
           continue
        end if
        ;= use permutation index of Gday_d
        ip_00 = Gsort_zero(ii_00,:,:); index of original day for the 99th percentile [nlat, nlon]
        ip = ip_00(j,i) ;day index
        iday= ind(ddd.eq.ip); indices of 24 time steps for Gday
        if(all(ismissing(iday))) then; skip some wrong indices
           continue
        end if
        Ta_Gday_zero(ii,:,j,i) = taC(iday,j,i)
        RH_Gday_zero(ii,:,j,i) = rh(iday,j,i)
        Rs_Gday_zero(ii,:,j,i) = Rs(iday,j,i)
        Rd_Gday_zero(ii,:,j,i) = Rd(iday,j,i)
        Rg_Gday_zero(ii,:,j,i) = Rg(iday,j,i)
        Ld_Gday_zero(ii,:,j,i) = Ld(iday,j,i)
        Lg_Gday_zero(ii,:,j,i) = Lg(iday,j,i)
        U_Gday_zero(ii,:,j,i)  = wind2m(iday,j,i)
        P_Gday_zero(ii,:,j,i)  = pa(iday,j,i)
        Z_Gday_zero(ii,:,j,i)  = zenith(iday,j,i)
        G_Gday_zero(ii,:,j,i) = Gzero(iday,j,i)


        ;=wet-bulb temperature
        TWday1d = TWday(ii_00,j,i)
        if(ismissing(TWday1d)) then; skip ocean grids
           continue
        end if
        ;= use permutation index of TWday
        ip_00 = TWsort_index(ii_00,:,:); index of original day for annual maximum  [nlat, nlon]
        ip = ip_00(j,i) ;day index 
        iday= ind(ddd.eq.ip); indices of 24 time steps for the day of TWday
        if(all(ismissing(iday))) then; skip some wrong indices
           continue
        end if
        Ta_TWday(ii,:,j,i) = taC(iday,j,i)
        RH_TWday(ii,:,j,i) = rh(iday,j,i)
        Rs_TWday(ii,:,j,i) = Rs(iday,j,i)
        Rd_TWday(ii,:,j,i) = Rd(iday,j,i)
        Rg_TWday(ii,:,j,i) = Rg(iday,j,i)
        Ld_TWday(ii,:,j,i) = Ld(iday,j,i)
        Lg_TWday(ii,:,j,i) = Lg(iday,j,i)
        U_TWday(ii,:,j,i)  = wind2m(iday,j,i)
        P_TWday(ii,:,j,i)  = pa(iday,j,i)
        Z_TWday(ii,:,j,i)  = zenith(iday,j,i)
        TW_TWday(ii,:,j,i) = TW3h(iday,j,i)

      end do
    end do

   end if

   ;================================
   ;---annual accumulative days with Gday >0 each year per grid cell

    Gday0 := where(Gday_sun.gt.0, 1, 0); daytime only (coszen>0)
    Gday0d_sun(ii,:,:) = tofloat(dim_sum_n(Gday0, 0))
    Gday0 := where(Gday_diff.gt.0, 1, 0);
    Gday0d_diff(ii,:,:) = tofloat(dim_sum_n(Gday0, 0))
    Gday0 := where(Gday_zero.gt.0, 1, 0);
    Gday0d_zero(ii,:,:) = tofloat(dim_sum_n(Gday0, 0))
    Gnight0 := where(Gnight.gt.0, 1, 0); night only (coszen<0)
    Gnight0d(ii,:,:) = tofloat(dim_sum_n(Gnight0, 0))
    TWday35 := where(TWday.gt.35, 1, 0); TWday>35
    TWday35d(ii,:,:) = tofloat(dim_sum_n(TWday35, 0))


  end do; end year loop


  ;====save output data=====
  ;define new coordinate variable
  time      :=  ispan(syr,fyr,1)
  time!0    = "time" ;named variables
  time@long_name = "simulation start:end years"
  time@units = "year CE"

  Gday_mx_sun!0="time"
  Gday_mx_sun&time = time
  Gday_mx_sun!1="lat"
  Gday_mx_sun&lat = lat
  Gday_mx_sun!2="lon"
  Gday_mx_sun&lon = lon

  copy_VarCoords(Gday_mx_sun, Gday_mx_diff)
  copy_VarCoords(Gday_mx_sun, Gday_mx_diff_L)
  copy_VarCoords(Gday_mx_sun, Gday_mx_diff_U)
  copy_VarCoords(Gday_mx_sun, Gday_mx_zero)
  copy_VarCoords(Gday_mx_sun, Gnight_mx)
  copy_VarCoords(Gday_mx_sun, TWday_mx)

  copy_VarCoords(Gday_mx_sun, Gday0d_sun)
  copy_VarCoords(Gday_mx_sun, Gday0d_diff)
  copy_VarCoords(Gday_mx_sun, Gday0d_zero)
  copy_VarCoords(Gday_mx_sun, Gnight0d)
  copy_VarCoords(Gday_mx_sun, TWday35d)


  ;define new coordinate variable for diel (8) time steps
  step      :=  ispan(1,nstep,1)
  step!0    = "step" ;named variables
  step@long_name = "time step of a day"
 
  Ta_Gday_sun!0="time"
  Ta_Gday_sun&time = time
  Ta_Gday_sun!1="step"
  Ta_Gday_sun&step = step
  Ta_Gday_sun!2="lat"
  Ta_Gday_sun&lat = lat 
  Ta_Gday_sun!3="lon"
  Ta_Gday_sun&lon = lon
 
  copy_VarCoords(Ta_Gday_sun, RH_Gday_sun)
  copy_VarCoords(Ta_Gday_sun, Rs_Gday_sun)
  copy_VarCoords(Ta_Gday_sun, Rd_Gday_sun)
  copy_VarCoords(Ta_Gday_sun, U_Gday_sun)
  copy_VarCoords(Ta_Gday_sun, Ld_Gday_sun)
  copy_VarCoords(Ta_Gday_sun, Lg_Gday_sun)
  copy_VarCoords(Ta_Gday_sun, P_Gday_sun) ; 
  copy_VarCoords(Ta_Gday_sun, Rg_Gday_sun)
  copy_VarCoords(Ta_Gday_sun, Z_Gday_sun) ; 
  ;copy_VarCoords(Ta_Gday_sun, Rsa_Gday_sun)
  copy_VarCoords(Ta_Gday_sun, G_Gday_sun) ; 

  copy_VarCoords(Ta_Gday_sun, Ta_Gday_diff)
  copy_VarCoords(Ta_Gday_sun, RH_Gday_diff)
  copy_VarCoords(Ta_Gday_sun, Rs_Gday_diff)
  copy_VarCoords(Ta_Gday_sun, Rd_Gday_diff)
  copy_VarCoords(Ta_Gday_sun, U_Gday_diff)
  copy_VarCoords(Ta_Gday_sun, Ld_Gday_diff)
  copy_VarCoords(Ta_Gday_sun, Lg_Gday_diff)
  copy_VarCoords(Ta_Gday_sun, P_Gday_diff) 
  copy_VarCoords(Ta_Gday_sun, Rg_Gday_diff)
  copy_VarCoords(Ta_Gday_sun, Z_Gday_diff) 
  copy_VarCoords(Ta_Gday_sun, G_Gday_diff) 

  copy_VarCoords(Ta_Gday_sun, Ta_Gday_zero)
  copy_VarCoords(Ta_Gday_sun, RH_Gday_zero)
  copy_VarCoords(Ta_Gday_sun, Rs_Gday_zero)
  copy_VarCoords(Ta_Gday_sun, Rd_Gday_zero)
  copy_VarCoords(Ta_Gday_sun, U_Gday_zero)
  copy_VarCoords(Ta_Gday_sun, Ld_Gday_zero)
  copy_VarCoords(Ta_Gday_sun, Lg_Gday_zero)
  copy_VarCoords(Ta_Gday_sun, P_Gday_zero)
  copy_VarCoords(Ta_Gday_sun, Rg_Gday_zero)
  copy_VarCoords(Ta_Gday_sun, Z_Gday_zero)
  copy_VarCoords(Ta_Gday_sun, G_Gday_zero)

  copy_VarCoords(Ta_Gday_sun, Ta_TWday)
  copy_VarCoords(Ta_Gday_sun, RH_TWday)
  copy_VarCoords(Ta_Gday_sun, Rs_TWday)
  copy_VarCoords(Ta_Gday_sun, Rd_TWday)
  copy_VarCoords(Ta_Gday_sun, U_TWday)
  copy_VarCoords(Ta_Gday_sun, Ld_TWday)
  copy_VarCoords(Ta_Gday_sun, Lg_TWday)
  copy_VarCoords(Ta_Gday_sun, P_TWday) 
  copy_VarCoords(Ta_Gday_sun, Rg_TWday)
  copy_VarCoords(Ta_Gday_sun, Z_TWday) 
  copy_VarCoords(Ta_Gday_sun, TW_TWday) 


 ;======save data as ncdf 
     dirNc  = dirPath3       ; where nc file will be saved
     filNc  = "G-daynight-"+syr+fyr+"-"+mo+"-vcbc-outdoor.nc"

     pthNc  = dirNc+filNc
     system("/bin/rm -f "+pthNc)       ; remove any pre-existing file
     ncdf = addfile(pthNc ,"c")        ; open output netCDF file

    ;===================================================================
    ; create global attributes of the file (optional)
    ;===================================================================
     fAtt               = True            ; assign file attributes
     fAtt@title         = "Body energy balance G "+syr+"-"+fyr

     fAtt@Conventions   = "None"
     fAtt@creation_date = systemfunc ("date")


     ncdf->Gday_mx_sun  = Gday_mx_sun
     ncdf->Gday_mx_diff = Gday_mx_diff
     ncdf->Gday_mx_diff_L = Gday_mx_diff_L
     ncdf->Gday_mx_diff_U = Gday_mx_diff_U
     ncdf->Gday_mx_zero = Gday_mx_zero
     ncdf->Gnight_mx  = Gnight_mx
     ncdf->TWday_mx  = TWday_mx

     ncdf->Gday0d_sun  = Gday0d_sun
     ncdf->Gday0d_diff = Gday0d_diff
     ncdf->Gday0d_zero = Gday0d_zero
     ncdf->Gnight0d = Gnight0d
     ncdf->TWday35d = TWday35d

     ;either define the time dimension UNLIMITED or make records with ncks
     system("echo o | ncks --mk_rec_dmn time "+pthNc+" "+pthNc)


     filNc0  = "Gday_G-Z-3h-"+syr+fyr+"-"+mo+"-vcbc-outdoor.nc"
     pthNc0  = dirNc+filNc0
     system("/bin/rm -f "+pthNc0)       ; remove any pre-existing file
     ncdf0 = addfile(pthNc0 ,"c")       ; open output netCDF file
     ncdf0->Z_Gday_sun = Z_Gday_sun ;
     ncdf0->Z_Gday_diff = Z_Gday_diff ;
     ncdf0->Z_Gday_zero = Z_Gday_zero ;
     ncdf0->G_Gday_sun = G_Gday_sun ;
     ncdf0->G_Gday_diff = G_Gday_diff ;
     ncdf0->G_Gday_zero = G_Gday_zero ;
     ncdf0->Z_TWday = Z_TWday ;
     ncdf0->TW_TWday = TW_TWday ;


     filNc1  = "Gday_sun-TQRUPZ-3h-"+syr+fyr+"-"+mo+"-vcbc-outdoor.nc" 
     pthNc1  = dirNc+filNc1
     system("/bin/rm -f "+pthNc1)       ; remove any pre-existing file
     ncdf1 = addfile(pthNc1 ,"c")       ; open output netCDF file

     ncdf1->Ta_Gday_sun = Ta_Gday_sun ; concurrent climate variables with Gday under sun scenario
     ncdf1->RH_Gday_sun = RH_Gday_sun
     ncdf1->Rs_Gday_sun = Rs_Gday_sun
     ncdf1->Rd_Gday_sun = Rd_Gday_sun
     ncdf1->Rg_Gday_sun = Rg_Gday_sun
     ncdf1->Ld_Gday_sun = Ld_Gday_sun
     ncdf1->Lg_Gday_sun = Lg_Gday_sun
     ncdf1->U2_Gday_sun = U_Gday_sun ; 2m wind
     ncdf1->P_Gday_sun = P_Gday_sun ;

     filNc2  = "Gday_diff-TQRUPZ-3h-"+syr+fyr+"-"+mo+"-vcbc-outdoor.nc"
     pthNc2  = dirNc+filNc2
     system("/bin/rm -f "+pthNc2)       ; remove any pre-existing file
     ncdf2 = addfile(pthNc2 ,"c")       ; open output netCDF file

     ncdf2->Ta_Gday_diff = Ta_Gday_diff ;concurrent climate variables with Gday under shade scenario
     ncdf2->RH_Gday_diff = RH_Gday_diff
     ncdf2->Rs_Gday_diff = Rs_Gday_diff
     ncdf2->Rd_Gday_diff = Rd_Gday_diff
     ncdf2->Rg_Gday_diff = Rg_Gday_diff
     ncdf2->Ld_Gday_diff = Ld_Gday_diff
     ncdf2->Lg_Gday_diff = Lg_Gday_diff
     ncdf2->U2_Gday_diff = U_Gday_diff 
     ncdf2->P_Gday_diff = P_Gday_diff 

     filNc4  = "Gday_zero-TQRUPZ-3h-"+syr+fyr+"-"+mo+"-vcbc-outdoor.nc"
     pthNc4  = dirNc+filNc4
     system("/bin/rm -f "+pthNc4)       ; remove any pre-existing file
     ncdf4 = addfile(pthNc4 ,"c")       ; open output netCDF file

     ncdf4->Ta_Gday_zero = Ta_Gday_zero ;concurrent climate variables with Gday under dark scenario
     ncdf4->RH_Gday_zero = RH_Gday_zero
     ncdf4->Rs_Gday_zero = Rs_Gday_zero
     ncdf4->Rd_Gday_zero = Rd_Gday_zero
     ncdf4->Rg_Gday_zero = Rg_Gday_zero
     ncdf4->Ld_Gday_zero = Ld_Gday_zero
     ncdf4->Lg_Gday_zero = Lg_Gday_zero
     ncdf4->U2_Gday_zero = U_Gday_zero
     ncdf4->P_Gday_zero = P_Gday_zero

     filNc3  = "TWday-TQRUPZ-3h-"+syr+fyr+"-"+mo+"-vcbc.nc"
     pthNc3  = dirNc+filNc3
     system("/bin/rm -f "+pthNc3)       ; remove any pre-existing file
     ncdf3 = addfile(pthNc3 ,"c")       ; open output netCDF file

     ncdf3->Ta_TWday = Ta_TWday ; simultanous climate variables at TWmax
     ncdf3->RH_TWday = RH_TWday
     ncdf3->Rs_TWday = Rs_TWday
     ncdf3->Rd_TWday = Rd_TWday
     ncdf3->Rg_TWday = Rg_TWday
     ncdf3->Ld_TWday = Ld_TWday
     ncdf3->Lg_TWday = Lg_TWday
     ncdf3->U2_TWday = U_TWday
     ncdf3->P_TWday = P_TWday

     ;either define the time dimension UNLIMITED or make records with ncks
     system("echo o | ncks --mk_rec_dmn time "+pthNc0+" "+pthNc0)
     system("echo o | ncks --mk_rec_dmn time "+pthNc1+" "+pthNc1)
     system("echo o | ncks --mk_rec_dmn time "+pthNc2+" "+pthNc2)
     system("echo o | ncks --mk_rec_dmn time "+pthNc3+" "+pthNc3)
     system("echo o | ncks --mk_rec_dmn time "+pthNc4+" "+pthNc4)


 end if





end 
